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1.0  INTRODUCTION 


The  purpose  of  this  report  is  to  evaluate  the  flutter  characteristics  of  a  control  fin  for  the 
Advanced  Kinetic  Energy  Missile  (AdKEM).  First  of  all  the  two-degree-of-freedom  equa¬ 
tions  of  motion  for  a  control  fin  under  aerodynamic  and  actuator  loads  will  be  determined. 

To  present  a  worst  case  analysis  of  the  problem,  these  equations  will  be  simplified  by  ne¬ 
glecting  mechanical  and  aerodynamic  damping  terms.  The  fin  characteristic  equation  will 
then  be  found  and  the  fin  flutter  regions  will  be  plotted.  Once  the  flutter  regions  are  plotted, 
the  fin  bending  stiffness  and  actuator  dynamic  stiffness  must  be  known  to  evaluate  the  possi¬ 
bility  for  fin  flutter  occurring. 

The  fin  bending  stiffness  is  determined  assuming  a  distributed  load,  equal  to  the  dynam¬ 
ic  pressure,  is  applied  to  the  fin.  The  fin  bending  moment  equation  is  integrated  numerically 
to  determine  the  fin  slope  and  resulting  bending  stiffness  at  any  point  along  the  fin  span. 
These  results  are  then  modified  for  input  into  the  stability  equation. 

The  flutter  analysis  contained  in  Appendix  A  assumes  that  the  actuator’s  hinge  moment 
stiffness  can  be  represented  by  a  second  order  system.  This  report  investigates  the  validity  of 
this  assumption  by  comparing  the  second  order  system’s  stiffness  to  that  of  a  third  order  sys¬ 
tem  and  a  detailed  non-linear  actuator  model.  Results  indicated  that  for  the  case  of  zero 
damping,  a  worst  case  flutter  condition,  the  second  order  system  model  provides  a  good  ap¬ 
proximation  of  the  real  actuator’s  dynamic  stiffness. 

The  values  found  for  fin  bending  stiffness  and  actuator  dynamic  stiffness  are  then  com¬ 
pared  to  the  fin  stability  plot  to  evaluate  the  possibility  of  fin  flutter.  The  results  show  that 
the  fin  design  considered  in  the  analysis  lies  in  the  flutter  region. 
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2.0  FIN  STABILITY  EQUATION 


The  equations  of  motion  of  a  missile  control  fin  may  be  determined  by  considering  the 
fin  free-body  diagram  as  shown  in  Figure  2-1.  This  figure  shows  a  two-degree-of-freedom 
control  fin  with  aerodynamic  and  actuator  loads  being  applied.  Appendix  A  contains  two 
reports  entitled  Two-Deeree-pf-Freedom  High  Speed  Flutter  Model  and  ARAP  Simplified 
Criterion  which  were  used  as  references  in  the  development  of  the  equations  of  motion. 
These  are 


Ixx  l\y 
IXy  lyy 


Q  +  (Q/Vm  A  CLatfyy  (Q/Vm  A  CU)fixy 
(Q/Vm  a  CL,.^xy  Qi  +  (Q/Vm  A  CLa)/3xx 


Q  A  CLo  Yac 

e 

Kd  +  Q  A  CU  Xac 

a 

e 

a 


(2-1) 


These  equations  can  be  simplified  by  assuming  that  all  damping  terms  are  zero.  This  results 
in  the  equations 


I XX 

Ixy 

Iyy 

Ke 

0 


Q  A  CL«  Yac  6 

IQ  +  Q  A  CLo  Xac  a 


(2-2) 


The  fin  characteristic  equation  can  now  be  determined  and  written  in  Laplace  form  as 


(Ixxlyy  Ixy  )s  +  (IxxKoa  +  lyyK#  —  IxyXafl)5^  +  KckjK#  =0  (2  —  3) 


where 

Kaa  =  Kd  +  Q  A  Cba  Xac  (2  —  4) 


Kafl  =  Q  A  CLo  Yac  (2-5) 

The  dynarr;c  pressure,  Q,  is  a  function  of  the  Mach  number  and  is  found  using 

0  -  j  na,,  P atm  Mach2  (2  —  6) 


The  lift  coefficient,  CLo  ,  center  of  pressure,  and  fin  surface  area  and  ineitias  of  the  fin 
shown  in  Figure  2-2  have  been  estimated  cs  shown  in  Table  2-1 
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TABLE  2-1.  Control  Fin  Characteristics 


MACH  NUMBER 

CLa(rad_1) 

Xac  (in) 

0.5 

4.0508 

-0.60 

1.0 

4.9217 

-0.78 

1.5 

4.5321 

0.15 

2.0 

3.3232 

0.28 

3.0 

2.1257 

0.26 

4.0 

1.5699 

0.204 

5.0 

1.2433 

0.17 

6.0 

1.0428 

0.15 

6.5 

0.9740 

0.145 

7.0 

0.8995 

0.14 

Fin  Reference  Area  «  3.463  in2 
Yac  =  0.9  in 

Ixx  =  285.43  (10-6)  in-lbr-s2 
Iyy  =  122.77  (10-6)  in-lbf-s2 
Ixy  =  28.672  (10-6)  in-lbf-s2 


The  location  of  the  center  of  pressure  in  the  x-axis,  X8C,  is  measured  from  the  fin  hinge 
line  and  is  positive  to  the  rear  as  shown  in  Figure  2-1.  The  roots  of  Equation  2-3  can  now 
be  solved  for  as  a  function  of  fin  bending  stiffness,  actuator  stiffness,  and  Mach  number  and 
the  stability  regions  of  the  control  fin  can  be  plotted.  A  listing  of  the  program  used  to  deter¬ 
mine  the  roots  of  Equation  2-3  is  provided  in  Appendix  B  and  the  results  are  plotted  in  Fig¬ 
ure  2-3. 
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3.0  CONTROL  FIN  BENDING  STIFFNESS 


The  fin  bending  stiffness,  Kq,  must  be  found  in  order  to  determine  if  the  control  fin  de¬ 
sign  is  unstable.  Bending  stiffness  is  defined  as 


K«(.V)  =  [dM-dfl]  |  y  (3-1 ) 

The  distance  y  is  measured  across  the  fin  span  from  the  root  as  indicated  in  Figure  2-1. 
The  equation  for  the  bending  moment,  assuming  a  constant  distributed  load,  Q,  applied  to  the 
fin,  may  be  found  by  integrating  the  shear  equation  for  the  steel  control  fin  root  shown  in 
Figure  3-1.  In  order  to  evaluate  the  shear,  the  fin  planar  area  must  be  determined.  This  is 

A(y)  =  J  y  [21,  +  (I,  - 1,)  y/ls]  (3-2) 


Now  the  shear  equation  can  be  found  by  considering  the  force  balance  at  a  point  y  along  the 
span. 


V(y)  =  Q[A(1S)-A(y)]  (3-3) 

where  A(ls)  is  the  total  fin  root  area.  Integration  of  Equation  3-3  results  in  the  bending  mo¬ 
ment  equation  as  a  function  of  y. 


M(y)  =  Q  [A(ls)  (y— Is)  1,  (y2  -  Is2)  -  0t  - 1,)  (y3-ls3)/61s]  (3— ») 


Next,  the  fin  slope  must  be  determined.  The  difference  in  slope  of  two  points  along  the 
fin  root  is  defined  as 

Yb 

eb  -  ea  =  /  (M/EI)  dy  (3-5) 

Y. 

At  the  fin  root  the  slope,  0a,  is  assumed  zero.  Therefore,  the  slope  at  any  location  y  may  be 
calculated  as 


Y 

0(y)  =  (1/E)  j  (M/I)dy  (3-6) 

o 


4 


assuming  that  the  modulus  of  elasticity,  E,  is  constant.  Since  both  the  bending  moment  and 
the  moment  of  inertia  are  functions  of  y.  Equation  3-6  could  be  better  solved  using  numeri¬ 
cal  integration.  In  that  case,  Equation  3-6  takes  the  form 

v 

(Ky)  =  (6y/E)  X  (3-7) 

j  =  o 


where  6y  is  the  integration  step  size.  In  order  to  determine  the  slope  at  a  point  on  the  fin,  the 
cross-sectional  moment  of  inertia  must  be  found  as  it  varies  with  span.  The  development  of 
these  equations  for  the  fin  root  in  Figure  3-2  is  shown  in  Appendix  C. 

Knowing  the  bending  moment  and  slope  at  any  location,  the  bending  stiffness  may  be 
calculated  from  Equation  3-1.  This  equation  can  also  be  expressed  as 


K0(y)  =  [M(y)  -  M(y-6y)]/[0(y)  -  6(y-6y)]  (3-8) 

where  the  notation  ‘(y-^y)’  denotes  the  moment  or  slope  at  the  previous  integration  point. 
Note  that  the  dynamic  pressure,  Q,  cancels  out  in  the  calculation  of  bending  stiffness. 

The  equations  discussed  above  have  been  incorporated  into  a  computer  program,  listed 
in  Appendix  D,  which  calculates  the  fin  bending  stiffness  and  other  parameters.  This  pro¬ 
gram  first  calculates  the  cross-sectional  moment  of  inertia  at  the  particular  span  location  us¬ 
ing  the  equations  developed  in  Appendix  C.  Then  the  bending  moment  at  that  location  is 
calculated  from  Equation  3-4.  This  is  used  in  Equation  3-7  to  calculate  the  slope  of  the 
point  and  the  bending  stiffness  is  determined  using  Equation  3—8.  Also  calculated  is  the  de¬ 
flection  of  the  fin  normalized  to  dynamic  pressure.  These  data  are  then  written  to  a  file  for 
plotting.  This  process  is  repeated  along  the  length  of  the  span. 

The  fin  geometry  parameters  shown  in  Figure  3-3  were  used  as  input  to  the  program 
discussed  above  and  the  results  are  shown  in  Figure  3-4.  This  plot  shows  the  fin  bending 
stiffness  and  normalized  deflection  as  they  vary  along  the  span. 

Since  the  fin  characteristic  equation  developed  in  Section  1.0  assumes  that  the  dynamic 
pressure  is  resolved  as  the  lift  force  applied  at  the  center  of  pressure,  a  comparison  of  the  two 
methods  has  to  be  made.  This  can  be  done  by  letting  the  deflection  of  the  distributed  load 
case  at  the  fin  tip  equal  the  deflection  of  the  simple  case  at  the  fin  tip  and  solving  for  the 
bending  stiffness  that  would  allow  this  deflection  to  occur.  This  is  shown  pictorially  in  Fig¬ 
ure  3-5.  Summing  moments  about  the  fin  root  for  the  simple  case  shows  that 


Fl  Y#c  =  Kes  6s  (3-9) 

Solving  for  the  spring  stiffness,  Kes,  by  relating  the  lift  force  to  dynamic  pressure  and  the 
bending  angle  to  deflection  shows  that 

Kes  =  [AOs)  Y,c  1s)/[z(ls)/Q]  (3-10) 
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where  z(ls)/Q  is  the  absolute  value  of  the  normalized  deflection  of  the  fin  at  the  tip  and  is 
taken  from  Figure  3—4.  Taking  the  value  of  z(ls)/0  from  this  figure  shows 


z(lsVO  =  -0.910(10-3') 


and  substituting  this  into  Equation  3-10  results  in  an  equivalent  spring  stiffness  of 


Kos  =  5280  in-lbf  rad 


This  value  can  now  be  used  in  conjunction  with  the  actuator  dynamic  stiffness  to  determine, 
using  Figure  2-3,  if  the  fin  design  lies  in  the  flutter  regions. 


4.0  ACTUATOR  DYNAMIC  STIFFNESS 

The  flutter  analysis  contained  in  Appendix  A  uses  a  second  order  linear  system  to  repre¬ 
sent  the  actuator’s  hinge  moment  stiffness.  The  intent  of  this  section  is  to  determine  how 
well  a  second  order  system  represents  the  dynamic  stiffness  of  an  actual  open  center  valve 
pneumatic  actuator.  To  accomplish  this  the  second  order  system’s  dynamic  stiffness  is  deter¬ 
mined  and  compared  to  that  of  a  linearized  actuator  model,  as  well  as  a  detailed  non-linear 
actuator  model. 

The  second  order  system’s  dynamic  stiffness  is  determined  by  considering  the  equation 
of  motion  of  a  second  order  system  subjected  to  a  torque  disturbance. 


Iyyri  +  Cad  +  K aa  =  T0cos(Qt) 

The  steady-state  solution  to  this  equation  is  given  by 

a  M  =  a*  cos(ftt-tf> )  (4-2) 

Substituting  Equation  4-2  and  its  derivatives  into  Equation  4-1  results  in 

-IyyGPa*  cos^t-^-QQa,  sin(£2t-0)  +  K act,  cos(Qt  -  <p)  =  Tocos(£2t)  (4-3) 

Since  each  term  in  this  equation  represents  a  force  acting  on  the  inertia,  the  equation  can  be 
represented  by  a  force  vector  polygon,  from  which  it  is  easily  seen  that 

T02  =  (Kofi,  -  IyyQ2at)2  +  (QjQ a,)2  (4-4) 

and 
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tan(0) 


C„Q/(Ka  -  lyyQ2) 


(4-5) 


Solving  for  the  amplitude  of  the  steady-state  fin  position,  aa,  and  recalling  expressions  for 
frequency  ratio,  natural  frequency  and  damping  factor  results  in 

«a  =  (T,,/K„)/[(l-r)2  +  {23rfp  (4-6) 


From  Equation  4-6  the  desired  expression  for  the  second  order  system's  dynamic  stiffness  is 
obtained 


Kd  =  OVaa)  =  K«  [(1-r2)2  +  (26r)2]!/2  (4-8) 


An  expression  for  the  dynamic  stiffness  of  the  third  order,  linearized  actuator  model, 
represented  by  the  block  diagram  in  Figure  4-1,  was  obtained  from  Reference  1  and  is  given 
by 


Kd  =  (1/s)  (IyyS3  +  Bs2  +  (Ka-Ha)  s  +  KaKaGc)  (4-9) 


Setting  s  equal  to  jQ,  an  expression  for  the  magnitude  of  Equation  4-9  is  obtained,  which  is 
the  desired  expression  for  the  third  order  system’s  dynamic  stiffness. 


Kj  =  {(K.  +  H0-I„q‘)2  +  [BQ-(K„K,Gc/Q)1j}* 


(4-10) 


Substituting  expressions  for  frequency  ratio,  natural  frequency  and  damping  factor,  which 
were  defined  for  the  second  order  system,  Equation  4-10  becomes 


Kd  =  K «  ([l-rMHa/Ke)]2  +  [^-(K.Gc/Q)])1 


(4-11) 


Note  that  for  the  condition  where  Ha  is  small  compared  to  Ka,  the  dynamic  stiffness  deter¬ 
mined  by  Equation  4-11  approaches  that  of  Equation  4-8  at  high  frequencies.  However,  at 
low  frequencies  the  dynamic  stiffness  determined  by  Equation  4-11  approaches  infinity, 
while  Equation  4-8  approaches  Ka  . 

In  order  to  apply  Equations  4-8  and  4-11,  an  expression  for  the  actuator’s  static  stiff¬ 
ness  must  be  determined.  From  the  diagram  shown  in  Figure  4-2,  an  expression  for  torque 
about  the  hinge  line  is  determined. 

T  =  [ACPC  -  (Ac  -  Ar)Ps]  L  (4-12) 
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Differentiating  this  equation  with  respect  to  a  results  in 


dT/da  =  AcL(dPc/da ) 


(4-13) 


From  the  actuator  geometry,  a  diftcrcntial  change  in  volume  is  related  to  a  differential 
change  in  actuator  position  accoiding  to  the  relation 

da/dV  =  -l/AfLscc-(a  )  (4-14) 


If  the  compression  of  the  gas  above  the  piston  follows  an  iscntropic  process  then  pressure 
and  volume  arc  related  according  to 


PcVn  =  constant 


(4-15) 


Differentiation  of  Equation  4-15  results  in 


dPc/dV  =  -nPc/V  (4-16) 

Combining  Equations  4-13,  4-14,  and  4-16  provides  an  expression  for  the  actuator’s  static 
stiffness. 


Ka  =  nPc(LAc)2sec2(a  )/V  (4-17) 

The  dynamic  stiffness  was  determined  for  the  small  and  large  AdKEM  pneumatic  ac¬ 
tuators  using  Equations  4-8,  4-11,  and  4-17.  The  following  parameters  were  used  for  the 
small  actuator. 

Pc  =  1587.3  psia  L  =  0.6  in 

V0  =  0.046  in3  n  =  1.667 

H,  =  0  in— Ibf/rad  Gc  =  1.0 

Iyy  =  122.77(10)-*  in-lbf-s2 

From  which  the  small  actuator's  static  stiffness  and  natural  frequency  are  determined  at  the 
null  position  (a  =  0°). 


Ac  =  0.2961  in2 
B  =  0.50  in-lbj/rad/s 
Kg  =  124.64  rad2/s 


K «  =  1808  in— lbf/rad 
=  611  Hz 


The  small  actuator’s  dynamic  stiffness  obtained  from  the  second  order  model  is  compared  to 
that  of  the  third  order  model  in  Figures  4-3  thru  4-5  for  various  values  of  6.  A  good  match 
is  seen  for  6  =  0.50,  which  requires  that  Q*  =  0.47.  A  comparison  for  the  case  of  zero 
damping  is  shown  in  Figure  4-6. 
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The  large  actuator's  dynamic  stiffness  was  determined  based  on  the  following  parame¬ 
ters. 

Pc  *  1492.2  psia  L  * 

V„  =  0  OS5  in ?  r,  = 

11,,  =  (I  in-lhf  rad  Cic  = 

Iyy  =  245.54(1(1)-*  in-lbrs- 

From  which  the  large  actuator's  static  stiffness  and  natural  frequency  arc  determined  at  the 
null  position  (a  =  ()*’) . 


0.560  in 
1  6<>7 
1.0 


Ac 

li 


0.586  in2 
0.50  in— lb|  rad 
124.04  rad2  s 


Ka  =  3138  in-lbf/rad 
=  569  Hz 


The  large  actuator's  dynamic  stiffness  obtained  from  the  second  order  model  is  compared  to 
that  of  the  third  order  model  in  Figures  4-7  thru  4-9  for  various  values  of  6.  A  good  match 
is  seen  for  6  =  0.25,  which  requires  that  Q,  =  0.44  .  A  comparison  fot  the  case  of  zero 
damping  is  shown  in  Figure  4-10. 

As  a  final  comparison  the  stiffness  of  the  small  AdKEM  actuator  was  determined  using 
a  detailed  non-linear  actuator  model.  A  block  diagram  of  this  model  is  shown  in  Figure 
4-11.  The  stiffness  was  determined  by  applying  a  torque  disturbance  and  determining  the 
resulting  deflection.  Simulation  results  are  shown  compared  to  the  second  and  third  order 
models  in  Figures  4-  4  and  4-6  for  6  =  0.50  and  6  =  0  respectively.  Note  that  for  the  case 
with  damping  the  detailed  non-linear  model  indicates  a  sharp  drop  in  stiffness  at  the  closed 
loop  system’s  natural  frequency,  which  is  not  predicted  by  the  linear  models.  However,  for 
the  case  of  zero  damping  the  models  are  in  fairly  close  agreement. 

Based  on  the  above  discussion  it  may  be  concluded  that  the  linear  second  order  system 
model  used  in  the  flutter  analysis  realistically  represents  the  stiffness  of  the  AdKEM  actuator 
for  the  case  of  zero  damping.  Further,  since  zero  damping  is  a  worst  case  for  flutter,  the  use 
of  a  second  order  system  model  with  zero  damping  should  provide  a  conservative  flutter  pre¬ 
diction. 
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5.0  CONCLUSIONS  AND  RECOMMENDATIONS 


Referring  to  the  values  of  control  fin  bending  stiffness  and  actuator  dynamic  stiffness 
determined  and  comparing  them  to  the  fin  stability  regions  shows  that  the  fin  lies  in  the  flut¬ 
ter  region  for  the  small  actuator  design.  This  is  shown  in  Figure  5  -1  It  is  also  seen  that  the 
large  actuator  design  is  very  close  to  the  flutter  regions.  One  way  of  eliminating  this  prob¬ 
lem  is  to  increase  or  decrease  the  actuator  stiffness.  Since  actuator  stiffness  is  largely  a  func 
tion  of  the  hinge  moment  requirement  this  value  can  not  be  readily  changed.  Therefore,  it  is 
required  that  the  fin  bending  stiffness  be  changed.  Referring  to  Figure  5-1  shows  that  the 
best  approach  would  he  to  increase  the  fin  bending  stiffness  because  this  would  avoid  the 
flutter  regions  should  the  hinge  moment  requiicmcnt,  and  thus  the  actuator  stiffness,  be  re¬ 
duced.  Therefore,  it  is  recommended  that  the  fin  root  thickness  be  increased  and  the  flutter 
analysis  be  performed  again. 
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Figure  2-1.  Two-Degree-of-Freedom  Control  Fin. 
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Figure  2-3.  AdKEM  Control  Fin  Stability  Regions. 
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Figure  3-4.  Control  Fin  Bending  Stiffness  and  Deflection  Curves 
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Figure  4-10.  AdKEM  Large  Actuator  Dynamic  Stiffness  for  6  =  0.0 
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NOMENCLATURE 


A  Control  fin  surface  area  fin2) 

A(y)  Control  fin  root  surface  area  as  it  varies  along  the  span  (ur) 
A<  Actuator  control  piston  area  (in2) 

Ar  Actuator  rod  area  (in2) 

B  Viscous  friction  coefficient  (in— lbf/rad/s) 

CLa  Control  fin  lift  coefficient  (rad~!) 

Ca  Actuator  rotational  damping  factor  (in-lbf-s/rad) 

Ce  Control  fin  bending  damping  factor  (in— Ibf — s/rad) 

E  Modulus  of  elasticity  for  steel  (lbf/in2) 

Fl  Lift  force  (lbf) 

Gc  Actuator  compensator  gain 

Ha  Hinge  moment  coefficient  (in— lbf/rad) 

Ixx  Fin  moment  of  inertia  about  x-axis  (in-lbf-s2) 

Ixy  Fin  product  of  inertia  (in-lbf-s2) 

Iyy  Fin  moment  of  inertia  about  y-axis  (in-lbf-s2) 

Ka  Equivalent  actuator  valve  gain  (rad/s/rad) 

K<j  Actuator  dynamic  stiffness  factor  (in-lbf/rad) 

Ka  Actuator  static  stiffness  (in-lbf/rad) 

Kq  Fin  bending  stiffness  factor  (in-lbf/rad) 

K0S  Simplified  fin  bending  stiffness  factor  (in-lbf/rad) 

L  Actuator  lever  arm  (in) 

lr  Fin  root  length  (in) 

ls  Fin  span  length  (in) 

lt  Fin  tip  length  (in) 

M(y)  Fin  bending  moment  as  it  varies  along  the  span  (in— lbf) 
Mach  Missile  Mach  number 

n  Specific  heat  ratio  of  helium 

najr  Specific  heat  ratio  of  air 

a  Actuator  position  (rad) 
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a  M  Actuator  steady-state  position  (rad) 

Palm  Atmospheric  pressure  (psia) 

Pc  Actuator  control  pressure  (psia) 

Ps  Actuator  supply  pressure  (psia) 

0  Dynamic  pressure  (psia) 

r  Frequency  ratio  (£2/Qn) 

pxx  Aerodynamic  viscous  damping  factor  (in-) 

Pxy  Aerodynamic  viscous  damping  factor  (in2) 

Pyy  Aerodynamic  viscous  damping  factor  (in2) 

T  Actuator  torque  output  (in-lbf) 

t  Time  (s) 

To  Torque  disturbance  (in-lbf) 

6  Damping  factor  (CaQn/2Ka) 

0  Fin  bending  angle  (rad) 

Q  Frequency  (rad/s) 

0(y)  Fin  slope  as  it  varies  along  the  span  (rad) 

Qn  Natural  frequency  (K^/Iyy)-?  (rad/s) 

6y  Numerical  integration  step  size  (in) 

V  Actuator  control  chamber  volume  (in3) 

Vm  Missile  velocity  (in/s) 

V(y)  Fin  shear  force  as  it  varies  along  the  span  (lbf) 
X8C  Location  of  center  of  pressure  from  x-axis  (in) 
<p  Phase  angle  (rad) 

y  Distance  along  fin  span  from  root  (in) 

Yac  Location  of  center  of  pressure  from  y-axis  (in) 
z(ls)  Fin  deflection  at  tip  (in) 
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APPENDIX  A 

TWO  DEGREE  OF  FREEDOM  HIGH  SPEED  FLUTTER  MODEL 
and  ARAP  SIMPLIFIED  CRITERION  REPORTS 


APPENDIX  A 

TWO  DEGREE  OF  FREEDOM  HIGH  SPEED  FLUTTER  MODEL 

Consider  a  fin  and  coordinate  system  located  at  fin  center  of  compliance. 


A  complete  summary  of  system  parameters  is  presented  in  Table  A-l. 

A  summary  of  the  dynamic  equation  of  motion  is  presented  in  Equation  A-l. 


A-l 


TABLE  A-l.  Description  of  System  Parameters 


Parameter 

Description 

Units 

A 

Fin  surface  area 

irr 

Ca 

Rotational  damping 

in-lb-sec/rd 

Ce 

Bending  damping 

in-lb-sec/rd 

CL  a 

Panel  lift  force  Coef. 

1/rd 

FP 

Lift  force  on  panel 

lbf 

Ixx 

Moment  of  inertia 
about  x  axis 

in-lb-sec2 

Iyy 

Moment  of  inertia 
about  y  axis 

in-lb-sec2 

Ixy 

Product  of  inertia 

in-lb-sec2 

K« 

Rotational  spring  rate 

in-lb/rd 

Ke 

Bending  spring  rate 

in-lb/rd 

M 

Bending  Moment 

in-lbf 

Q 

Dynamic  Pressure 
(l/2pV2) 

lb/in2 

T 

Shaft  Torque. 

in-lb 

V 

Velocity 

in/sec 

Xac 

Dist.  of  C.P.  from 
x  axis 

in 

Xcg 

Dist.  of  C.G.  from 
x  axis 

in 

Yac 

Dist.  of  C.P.  from 
y  axis 

in 

Ycg 

Dist.  of  C.G.  from 
y  axis 

in 

a 

Rotational  angle 

rd 

Y 

Ixy/  yAxx  Iyy 

1 

e 

Bending  angle 

rd 

rixx 

Aero  viscous  damping 
Cocff. 

in2 

*lxy 

Aero  viscous 
damping  Coeff. 

in2 

%y 

Aero  viscous 
damping  Coeff. 

in2 

Uy 

Ixv  lyy 


a 

L  J 


4  W  ‘  [?  4  [-']  •  [“] 


(A-l ) 


The  aerodvnamic  moments  presented  i’i  Equation  A-l  may  be  expressed  b\ 
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Equation's  A-l  and  A-2  may  be  combined  to  vield: 
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Simplifying  Equation  A-3  Produces: 


V  A  CLa)  ^ 

[:  =]}  M 


*  0 


(A-3) 


Q,  +  (Q/V  A  CLa)  fjyy 

(Q/V  A  CLa)  Vxy  e 

(Q/V  A  CLa)  rjr, 

Ca  +  (Q/V  A  CU)  [« 

Kfl 

0 


Q  A  CLaYac 
(Kj  +  Q  A  CLaXac) 


0  (A— 4) 


Consider  solving  for  the  roots  of  the  characteristic  equation  of  expression  4  by  taking 
the  Laplace  transform  of  Equation  A-l  This  expression  may  then  be  expressed  as: 


(Ixx  S2  +  CgS  +  Kg)  (IXy  S2  +  Cg^S  +  Kafl)  1  \e] 

=  0  (A-5) 

Where: 

(I*y  S2  +  C^j)  (lyy  S2  +  Car  s  +  Kaa)J 

=  Q,  +  (Q/V  A  CLa)  TJyy 

(A-6a) 

=  (Q/V  A  CLa)  >7xy 

(A-6b) 

Cd 

=  Ox  +  (Q/V  A  CLa)  >7 xx 

(A-6c) 

A-3 


0  A  CLqYhc 


(A~6d) 


K a*  * 


Kaa  =  Ka  +  0  A  CUXnc  lA_(>c) 

The  systems  characteristic  equation  is: 

(Ixx  ::2  Cfl  s  +  K$)  (Ivy  s*  +  Q,  s  +  K,,,,)  - 
(Ixy  s2  +  Ca0  s  +  V^o)  (Ixy  s2  +  Ca0  s)  =  0 

Which  is  also: 

(Ixx  Iyy  -  Ixy2)  s4  +  (ixx  Ca  +  Iyy  C0  -  2Ixy  Ca0 )  s3 
+  (Ixx  Kaa  +  C0Ca  +  Iyy  Ko  ~  Ca02  -  Ixy  K^)  s2 
+  (Kaa  Cj  +  K</Q  —  K aO  S  K#  Kaa  ~  0 

Equation  A- 7  is  the  characteristic  equation  from  which  the  appropriate  roots  may  be 
computed  to  determine  stability. 

The  technique  for  studying  flutter  consists  of  finding  the  flight  conditions  which  gener¬ 
ate  roots  of  Equation  A-7  with  positive  real  parts.  These  solutions  generate  divergent  oscil¬ 
lations.  A  sample  solution  for  a  high  v  flutter  problem  is  found  in  Figure  A— 2.  This  figure 
shows  the  boundary  between  flutter  and  flutter  free  operation  as  a  function  of  rotational 
springrate,  Ka  ,  and  rotational  damping,  Co  .  Note  that  Ka  and  Q,  are  dominated  by  actua¬ 
tor  springrate  and  damping.  Note  also  that  the  flutter  boundary  is  plotted  for  five  values  of 
bending  springrate  at  a  constant  Mach  number  and  altitude.  Figure  A-2  illustrates  that  there 
are  three  ways  to  prevent  simple  fin  flutter. 

1.  Generate  a  high  value  of  rotational  springrate,  Ka .  Fcr  Ka  >  8200  in-lb/rd  no  combi¬ 
nation  of  conditions  can  generate  flutter. 

2.  Generate  a  high  value  of  bending  springrate,  K e  ■  For  K$  >  10000  in-lb/rd  flutter  can 
not  occur.  There  is  no  value  of  Ka  and  Q,  that  will  generate  flutter. 

3.  Generate  sufficient  K#  and  Ca  to  prevent  actuator  stiffness  and  damping  characteristics 
from  intersecting  the  flutter  zone.  Note  a  hydraulic  damper  may  be  required  to  achieve 
sufficient  rotational  damping  in  such  a  scheme. 

Note  that  the  technique  number  one  is  the  traditional  solution  to  fin  flutter.  However, 
techniques  two  and  three  may  be  required  with  a  pneumanic  fin  actuation  system. 

Figures  A-3  and  A-4  have  been  included  to  show  the  affect  of  damping  on  the  flutter 
boundaries.  The  analysis  of  Figure  A-3  shows  where  the  flutter  boundaries  exist  when  the 
bending  damping,  Q  ,  is  reduced  to  0.  Likewise,  Figure  A-4  shows  the  affect  of  reducing 
the  aerodynamic  damping  terms  to  50%  of  the  value  calculated  for  the  specified  flight 
condition. 
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A  RAP  SIMPLIFIED  CRITERION 


A  simplified  criterion  was  developed  by  GPSD  and  Aeronautical  Research  Associates 
of  Princeton  in  order  to  evaluate  IQ,  and  IQ?  for  flutter  free  operation.  This  criterion  only 
applies  to  high  v  supersonic  t  ins  A  summary  of  the  criterion  follows. 

If  all  damning  terms  are  neglected  in  Equation  A—!  then: 


lxx  Ixy 
Ixy  lyy 


Ko  Q  A  CL*,  Vac 
0  (IQ,  +  0  A  CLa  Xac) 


Likewise  the  system  characteristic  equation  becomes: 


0 

(A— 8) 


(lxx  lyy -Ixy7)  s4  +  (ix\  Koa  +  tyy  K®-Ixy  Kao)  s2  +  Kaa  IQ, 


0  (A-9) 


If  the  characteristic  equation  is  factored,  the  condition  where  all  four  roots  have  the 
same  frequency  (that  is  the  condition  where  flutter  begins)  is: 


(lxx  Koa)  -  2  lxx  lyy  K aa  *0?  -  2  lxx  Ixy  IQ,a  K a0 

-  2  lyy  Ixy  IQ,  IQ,*  +  (lyy  Ko)2  +  (Ixy  Ka0)2  (A-10) 

+  4  Ixy2  IQ,  Koa  =0 

Equation  A-10  which  is  the  boundary  of  flutter  is  plotted  for  the  sample  case  in  Figure 
A-5.  In  this  figure  the  fin  flutter  boundary  is  presented  as  a  function  of  rotational  and  bend¬ 
ing  stiffness.  Of  particular  interest  are  points  A  and  B  on  this  curve.  If  bending  stiffness, 

IQ, ,  is  greater  than  that  defined  by  point  A  then  no  flutter  can  occur.  If  rotational  stiffness, 

K a ,  is  greater  than  that  defined  by  point  B,  men  no  flutter  can  occur.  Points  A  and  B  are  de¬ 
fined  by: 

Point  A  IQ,  >  Ixx/Ixy  Q  A  CLq  Yac 

(No  Flutter)  (A-lla) 

Point  B  Ka  >  Iyy/Ixy  Q  A  CLa  Xac  -  Q  A  Cl*,  Xac 

(No  Flutter)  (A-llb) 

Refering  to  Figures  A-2,  A-3,  A-4,  and  A~5,  the  following  conclusions  may  be 
reached: 

1.  Figure  A-2  illustrates  that  increasing  Q,  can  increase  the  value  of  IQ,  that  is  required 
for  absolute  stability.  Figures  A-3  and  A-4  illustrates  that  decreasing  either  bending 
damping  or  aerodynamic  damping  will  increase  the  rotational  stiffness  necessary  for 
absolute  stability  (Point  B). 

2.  The  analysis  of  Figures  A-2,  A-3,  and  A-4  shows  that  the  criterion  for  the  bending 
springrate  required  for  absolute  stability,  (Point  A),  is  affected  only  slightly  by  system 
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damping.  The  bending  stiffness  requirement,  for  absolute  stability,  that  is  predicted  as 
suming  zero  damping  is  9650  in-lb/rd.  Variations  in  the  damping  terms  increased  this 
figure  only  slightly  to  10,000  in-lb/rd. 

The  simple  criterion  described  above  can  therefore  be  used  as  a  useful  guide  to  predict 
the  bending  stiffness,  K.u.  that  can  be  used  to  eliminate  flutter  regardless  of  the  values  of  ro¬ 
tational  stiffness,  Kq  .  and  rotational  damping,  Ca  . 

A  more  detailed  analysis  is  required  to  predict  the  rotational  stiffness,  ,  that  is  re¬ 
quired  to  produce  absolute  stability  because  of  the  significant  aflects  of  system  damping  on 
the  flutter  boundaries.  Using  such  an  approach  to  eliminate  flutter  is  therefore  much  more 
involved  and  may  requite  controlling  both  rotational  stiffness,  K<, ,  and  damping,  Ca  . 
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Bending  Springrate  I  xx  =  3.58  lb  in2 

0  KTH  =  5000  I  yy  =  2.92  lb  in2 

o  KTH  =  6000  I  xy  =  1.64  lb  in2 
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Figure  A-2. 


I  xx  =  3,58  lb  in2 
I  yy  =  2.92  lb  in2 
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STABILITY  EQUATION  PROGRAM  LISTING 

PROGRAM  NAME:  FINS1MP  FOR 


Steve  Cayson 

US  Army  Missile  Command 
Research,  Development,  and  Engineering  Center 
Redstone  Arsenal,  AL 
OCT-NOV  89 


This  program  is  used  to  calculate  the  flutter  regions  for  the 
AdKEM  control  section  fin  using  equations  developed  in  a  paper 
written  by  Garrett  Fluid  Systems  Division  in  Phoenix,  A Z.  The 
equations  were  developed  using  a  two-degree-of-freedom  model  for 
a  control  surface  assuming  no  damping  terms.  The  characteristic 
equation  for  this  system  was  then  found  and  used  to  calculate  the 
fin  stabilty  regions  based  on  fin  and  actuator  stiffness  and 
aerodynamic  parameters.  The  region  is  solved  for  and  the  data  is 
written  to  a  file,  FLUTTER.DAT,  that  can  be  plotted. 


VARIABLE  DESCRIPTION 


IXX 

IYY 

IXY 

REFA 

Q 

MACH 

V 

CLIFT 

XCP 

YCP 

KALPHA 


Fin  inertia  about  hinge  line  (in-lbf-s2) 

Fin  inertia  about  root  line  (in-Ibf-s2) 

Fin  product  of  inertia  (in-lbf-s2) 

Control  fin  surface  area  (in2) 

Dynamic  pressure  (psia) 

Mach  number 
Missile  velocity  (in) 

Fin  lift  coefficient  (rad-1) 

Distance  of  center  of  pressure  from  x-axis  (in) 
Distance  of  center  of  pressure  from  y-axis  (in) 
Fin  bending  stiffness  factor  (in-lbf/rad) 


REAL  IXX,  IYY,  IXY,  KAA,  KAT,  KALPHA,  KTHETA,  MACH 


DIMENSION  A(ll),  B(ll),  C(26),  D(26),  E(25),  F(25) 

DIMENSION  RR(1G),  RI(iO)  ,AA(26),  BB(26),  RREAL(25),  RIMAG(25) 
OPEN(4,  FILE=’FLU1TER.DAT’  ,STATUS='NEW’) 

*****  CONSTANTS  USED  IN  ROOT  SOLVER  SUBROUTINE  ***** 


-  ORDER  OF  POLYNOMIAL  EQUATION 
NN  =  2 
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C  -  INITIAL  ESTIMATES  FOR  ROOTS 
XO  =5. 

YO  =  5. 

X  =  XO 

V  =  YO 

C . INPUT  FIN'  CHARACTERISTICS . 

C  -  INERTIAS  (in-lbf-s2) 

IXX  =  .00028543 
IYY  =  .00012277 
IXY  =  .000028672 

C  -  FIN  REFERENCE  AREA  (in2) 

REFA  =  3.463 

C  -  INPUT  MACH  NUMBER  AND  CALCULATE  DYNAMIC  PRESSURE  (psia) 
WRITE(*,  100) 

100  FORMAT(2X,  ’MACH  NUMBER?’, 2X) 

READ  (*,  110)MACH 
110  FORMAT(  F10.7) 

Q  =  10.29*  (MACH* *2.) 

C  -  CALCULATE  LOCAL  VELOCITY  (in/s)  ASSUMING  TEMP  70F 

V  «  MACH*  13540. 

C  -  INPUT  FIN  LIFT  COEFFICIENT  (rad"1)  AND  CP  (in)  AT  MACH  NUMBER 
WRITE(*,  120)MACH 

120  FORMAT(/,2X,’LIFT  COEFFICIENT  AT  MACH  \  F5.2,’  ?  (l/rad)\  2X) 
READ(*,  110)  CLIFT 
WRITE(*,  130) 

130  FORMATS  2X,  ’DIST  OF  CP  FROM  HINGE  LINE  (+  is  in  aft  direction)’ 

&,  2X) 

READ(*,  110)XCP 
YCP  =  0.91 

C  *****  LOOP  TO  DETERMINE  FLUTTER  REGION  ***** 

DO  200  K  =  0, 10000,  20 
KALPHA  =  K*1.0 

C  -  PRELIMINARY  CALCULATIONS 
QUANT  =  Q*  REFA*  CLIFT 
KAT  =  QUANT*YCP 
KAA  =  KALPHA  +  QUANT*XCP 

C  -  REAL  COEFFICIENTS  OF  POLYNOMIAL  -  A(l)  IS  HIGHEST  ORDER  COEF. 
A(l)  =  IYY**2. 

A(2)  =  4*KAA*(IXY**2.)  -  2*IYY*IXY*KAT  -  2*IXX*IYY*KAA 

A(3)  =  (IXX*IXX*KAA*KAA)  +  (IXY*KAT)**2.  -  2* IXX* IXY* KAA* KAT 
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C  -  IMAGINARY  COEFFICIENTS  OF  POLYNOMIAL 
B(l) *  0. 

B(2)  =  0. 

B(3) =  0. 

C  ~  SOLVE  FOR  ROOTS  USING  SUBROUTINE  ROOTS 
CALL  ROOTS(A,  B,  NN,  X.  Y,  RR.  RI,  NERR) 

C  -  CHECK  TO  SEE  IF  ROOTS  ARE  EQUAL,  IF  SO  HALT  PROGRAM  EXECUTION 
CHECK  =  RR(1)  -  RR(2) 

IF  (ABS(CHECK)  .LT.  .1)  GOTO  99 

IF  (RR(1)  .LT.  0.0)  RR(1)  =  0.0 
IF  (RR(2)  .LT.  0.0)  RR(2)  =  0.0 

C  WRITE(MO)  KALPHA,  RR(1),  RR(2) 

WRITE(4,  10)  KALPHA,  RR(1),RR(2) 

10  FORMAT(2X,  F7.I,  2(2X,  FI  0.3)) 

200  CONTINUE 
99  END 


SUBROUTINE  ROOTS(A,  B,  NN,  X0,  Y0,  RR,  RI,  NERR) 

DOUBLE  PRECISION  AA,  BB,  RREAL,  RIMAG,  C,  D,  E,  F,  P,  Q,  X,  Y,  TOL,  DENOM, 
R  RS 

DIMENSION  A(*),  B(*),  RR(*),  RI(*) 

DIMENSION  AA(26\  BB(26),  C(26),  D(26) 

DIMENSION  RREAL(25),  RIMAG(25),  E(25),  F(25) 

C  INITIALIZE 
NERR  =  0 
R  =  0.0D0 
N=  NN 
NP1  =  N+l 
K  =  0 

TOL  =  .5D-8 
X  =  X0 
Y  =  Y0 

DO  5  1=1,  NP1 
AA(I)  =  A(I) 

5  BB(I)  =  B(I) 

C  TEST  FOR  FIRST  DEGREE  EQUATION 
10  IF(N  .EQ.  1)  GOTO  60 

C  BEGIN  SYNTHETIC  DIVISION  FOR  F(Z)  AND  F’(Z) 

20  KTR  =  0 
NP1  =  N+l 

IF(X  .EQ.  0.0D0)  X=.37D0 
IF(Y  .EQ.  0.0D0)  Y=.37D0 
C(l)  =  AA(1) 

D(l)  =  BB(1) 

E(l)  =  C(l) 
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F(l)  *  D(l) 

23  DO  25  I=2,NP1 
TM1  =  1-1 

C(I)  =  AA(1)  +  X*C(IMI)  -  Y*D(IM  1 ) 

I)(I)  =  BIJ(I)  +  Y*C(IM  1 )  +  X*l)(IM ! ) 

IF(I  .EO.  NP1)  GOTO  25 

E(I)  =  C(l)  +  X*E(IM1)  -  Y * F( I M 1 ) 

F(l)  =  D(I)  -i-Y*E(IMl)  +X*F(1M1) 

25  CONTINUE 

C  CALCULATE  X  AND  IY 

30  DENOM  =  E(N)**2+F(N)**2 

P=  (C(NP1)*E(N)  +  D(NPl)*F(N))/DENOM 
Q  =  (D(NP1)*E(N)  -  C(NPl)*F(N))/DENOM 
X  =  X-P 

Y  =  Y-Q 
RS  =  R 

R  =  DSQRT(P**2+Q**2) 

IF(R  .LT.  RS)  GOTO  40 

C  TEST  LOOP  COUNTER  WHEN  CORRECTION  DOES  NOT  DECREASE 

KTR  =  KTR+1 
IF(KTR  .LE.  25)  GOTO  40 
NERR  =  1  GOTO  50 

C  CHECK  FOR  CONVERGENCE 

40  IF(R.GT.TOL)  GOTO  23 
C  ROOT  FOUND  -  REDUCE  POLYNOMIAL 

50  K  =  K+l 

RREAL(K)  =  X 
RIMAG(K)  =  Y 
DO  55  1=1,  N 
AA(I)  =  C(I) 

55  BB(I)  =  D(I) 

N  =  N-l 
GOTO  10 

C  SOLUTION  FOR  FIRST  DEGREE  POLYNOMIAL 

60  DENOM  =  AA(1)**2+BB(1)**2 

x  =  (-AA(1)*AA(2)  -  BB(l)*BB(2))/DENOM 

Y  =  (  AA(2)*BB(1)  ~  AA(l)*BB(2))/DENOM 
K  =  K+l 

RREAL(K)  =  X 
RIMAG(K)  =  Y 

C  MOVE  DOUBLE  PRECISION  ROOTS  TO  SINGL  E  PRECISION  OUTPUT 

70  DO  75  1=1, NN 
RR(I)  =  RREAL(I) 

75  RI(I)  =  RIMAG(I) 

C 

999  RETURN 
END 
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APPENDIX  C 

FIN  MOMENT  OF  INERTIA  EQUATIONS  DEVELOPMENT 


APPENDIX  C 


The  leading  wedge  of  the  control  fin  root  is  represented  as  shown  in  Figure  3-2.  The 
width.  h3(y).  is  calculated  using 


My;  =  (lo  -  If* Ky  is)  +  lt.i  (C-l  ) 

and  the  thickness,  t3(v).  is  found  from  the  equation 

t.i(y)  =  (ts-tf)(y/ls)  +  tr  (C-2) 

The  leading  edge  thickness  can  be  calculated  using 

tlc(y)  =  Oslc  —  hlc)(y/^s)  hie  (C— 3) 

A  trapezoid  such  as  the  one  shown  below  has  now  been  defined  and  the  inertia  about 
the  x-axis  can  be  found. 


e(Y) 

i 

_ i 

r~ - 

-  t3( 

Figure  C-l.  Typical  Fin  Wedge  Cross-Section 

The  trapezoid  may  be  divided  into  two  regions,  the  rectangular  portion  and  the  triangular 
portion.  The  inertia  of  the  rectangular  portion  is  just 


l3R(y)  =  [bj(y)  n«(y)3]/l2  (C-4) 

The  inertia  of  the  triangular  sections  can  be  found  by  integration  of 

I3T  =  2J  z2  dA  (C-5) 

A 

where  z  is  the  distance  as  shown  in  Figure  C-l.  Substitution  of  the  proper  limits  into  Equa¬ 
tion  C-5  results  in  the  double  integral 

bj(Y)  mux+bu 

I3T  =  2  j  j  z2  dzdx  (C-6) 

0  ba 


-1 


where 


bz.i  =  i|e(v  )/2 


an  J 


(C-7) 


=  |t3(y )  —  t|f(y)|  '(2  b3(y))  (C-8) 

Evaluation  of  Equation  C-fc  tesults  in 

I.rr(>')  =  bi(y)2  {mz.rV,(y)2/6  +  2mZ32b/3b3(y)/3  +  m^b^2]  (C-9) 

and  the  total  inertia  of  the  ieading  wedge  section  as  a  function  of  y  is 

h(y)  =  i3R(y)  +  I3r(y)  (C-iO) 

The  inertia  of  the  trailing  section  is  found  similarly.  The  inertia  of  the  rectangular  por¬ 
tion  is 

im(y)  =  (bi(y)  iie(y)3]/i2  (c-ii) 

and  the  inertial  of  the  triangular  section  is 

Ii-r(y)  =  [bi(y)2  [mzl3b]fy)2/6  +  2m*i2bzibi(y)/3  +  mzlbzi2]  (C— 12) 

where 


bi(y)  =  (1,1  -lri)(y/ls)  +  ln  (C— 13) 

ti(y)  =  (ts  -tr)(y/ls)  +  tr  (C-14) 

he(y)  =  Oste  —  Irte)  (y/ls)  +  Irte  (C— 15) 

bzt  =  tle(y)/2  (C— 16) 

mzt  =  [ti(y)-tle(y)]/[2b1(y)]  (C-17) 

Therefore,  the  total  inertia  of  the  trailing  wedge  is 

ii(y)  =  iiR(y)  +  My)  (c-18) 

The  middle  section  of  the  fin  is  simply  a  rectangle  and  the  inertia  is  found  using 

I2(y)  =  b2(y)  t2(y)3/l2  (C-19) 

C-2 


where 


b;(y)  =  ( 1 12  -  lr:)  (>'/ls^  +  lr2  (C-20) 

and 

t’(y)  =  (t*  — »,)  (y/is)  +  tr  (C-2 1 ) 

The  total  inertia  of  the  section  is  then  found  from 

1.0, (V)  =  Ii(y)  +  I  _■>(>’)  +  l3(y)  (C-22) 


C-3/(C-4  Blank) 
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BENDING  PROGRAM  LISTING 


PROGRAM  NAME  Fl’LLFIN  FOR 


Steve  Cayson 

l '5  Army  Missile  Command 

Research,  Development,  and  Engineering  Center 

Control  Systems  Branch 


04  DEC  89 


This  program  will  calculate  the  fin  bending  spring  rate  for  the 
stainless  steel  core  of  the  AdKEM  control  fin  assuming  a  constant 
distributed  load  (dynamic  pressure)  on  the  fin.  First  of  all  the 
fin  geometry  is  input  to  the  program.  This  includes  the  span, 
various  thicknesses,  and  the  root  and  tip  lengths.  Then  a 
loop  is  used  to  calculate  the  moment  and  slope  at  each  point  along 
the  span  and  solve  for  the  stiffness  measured  in  in— lbf/rad. 

The  equation  is 


K(Y)  =  dMOMENT/dTHETA 
VARIABLE  DESCRIPTION 


THKFACT 

LR1 

LR2 

LR3 

LRTOT 

TROOT 

LT1 

LT2 

LT3 

LTTOT 

TSPAN 

LSPAN 

TSLE 

TRLE 

TSTE 

TRTE 

ATOT 

MOM1 


Thickness  multiplier  for  variable  thicknesses 
Trailing  wedge  length  at  root  (in) 
Mid-section  length  at  root  (in) 

Leading  wedge  length  at  root  (in) 

Total  fin  length  at  root  (in) 

Mid-section  thickness  at  root  (in) 

Trailing  wedge  length  at  tip  (in) 

Mid-section  length  at  tip  (in) 

Leading  wedge  length  at  tip  (in) 

Total  fin  tip  length  (in) 

Mid-section  thickness  at  span  (in) 

Length  of  span  (in) 

Thickness  of  leading  edge  at  span  (in) 
Thickness  of  leading  at  root  (in) 

Thickness  of  trailing  edge  at  span  (in) 
Thickness  of  trailing  edge  at  root  (in) 

Total  fin  root  area  (in2) 

Bending  moment  at  fin  root  (in-lbf) 


IMPLICIT  REAL  (I  ,K-M) 

DOUBLE  PRECISION  THETA1,  THETA2,  DELTHT,  XDEF,  DY 
OPEN(4,  FILE=’FULL.DAT’,  STATUS  =  ‘NEW’) 

THKFACT  =  1.0 


. INPUT  FIN  GEOMETRY  SPECIFICATIONS  ***** 
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o  n 


C  ROOT  GEOMETRY  (ALL  DIMENSIONS  IN  INCI IES) 
LR1  =  .699 
LR2  =  .822 
LR3  ^  .674 

LRTOT  =  LRI  ^  LR2  -  1.R3 
TROOT  =  THKFACT  •  0.1 


C  TIP  GEOMETRY 
LT1  =  .299 
LT2  =  .444 
LT3  =  .357 

LTTOT  =  LTI  +  LT2  +  LT3 
TSPAN  =  THKFACT  *  0.057 

C  SPAN  LENGTH 

LSPAN  =  1.80 


C  LEADING  AND  TRAILING  EDGE  THICKNESSES  AT  ROOT  AND  SPAN 
TSLE  =  THKFACT  *  .017 
TRLE  =  THKFACT  *  .025 
TS'TE  =  THKFACT  *  .025 
TRTE  =  THKFACT  *  .025 

C  CALCULATE  TOTAL  FIN  AREA  AND  BENDING  MOMENT  AT  Y=0 
ATOT  =  LSPAN/2.  *  (LRTOT  +  J.TTOT) 

MOM1  =  LSPAN/6.  *  (LSPAN*(LTTOT  +  2.*LRTOT)  -  6.*ATOT) 

C  . OTHER  VARIABLES  ***** 

C  FIN  COEFFICIENT  OF  ELASTICITY  (psi) 

E  ~  28.5E06 

C  INTEGRATION  STEF  SIZE 
DY  =  .001 

C  DO  LOOP  TERMINATION  VALUE  CALCULATION 
JOY  =  LSPAN/DY 

C  LOOP  INITIALIZATION 
SUMMOI  =  .0 
THETA  1  =  .0 
ZDEF  =  .0 
JPRINT  ^  1 


C 

C 

c 


!!!!!!!!!!!!!!!!!! !!!!!!!!!!!!! i !!!! H ! 1 < 1 !!!!!!! i  J ! !!!!!!!!!!!!!!! !!!!!!!!i !!!!!!!!!!!!!!!!!!!!! !!!!!!!!!! 
!!!.'!!!!!!!!!!!!!!!  [DON’T  MESS  WITH  ANYTHING  BFLOW  HERE!!:!!!!!!!!!!!!!!!! 
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!S!!!!!!!!!!!!!!!!!!!!!!!!!j!!S!!!! !!!!!!!!!:!!!!!! !!!!!!! 
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C  *****  BEGIN  NUMERICAL  INTEGRATION  TO  SOLVE  FOR  STIFFNESS  ***** 
DO  10  J  =  1,  JDY-9 
Y  =  J*DY 


*****  CALCULATE  CROSS-SECTIONAL  INERTIA  OF  FIN  AT  Y  *  *  *  *  * 
CALCULATE  TfLVILING  SECTION  PARAMETERS 
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B1  =  (LT1~LR1)*Y/LSPAN  +  LR1 
T1  =  (TS  PAN -TROOT)  *Y/LSPAN  +  TROOT 
TTE  =  (TSTE-TRTF)  *  Y/LSPAN  +  TRTE 
MZ1  =  (Tl-TTE  )/(2.*Bl ) 

BZ1  =  TTE '2. 

C  CALCULATE  TRAILING  SECTION  INERTIA 
C  (RECTANGULAR  PORTION) 

XII R  =  (1/12. )*B1  *(TTE**3) 

C  (TRIANGULAR  PORTION) 

XI1T  =  Bl*(((MZl*Bl)**3)/6.  +  (2.*BZl/3.  )*(MZ1  *B1)*  *2 + 

&  MZ1*B1*BZ1**2) 

C  TOTAL 

XU  =  XUR  +  XIIT 

C  CALCULATE  MIDDLE  SECTION  INERTIA 
Cl  =  (LT2-LR2)  *  Y /LS  PAN  +  LR2 
C2  ”  ((TSPAN-TROOT)*  Y/LSPAN  +  TROOT)*  *3 
C3  =  1/12.  XI2  =  C3*C1*C2 

C  CALCULATE  LEADING  SECTION  PARAMETERS 
B3  =  (LT3-LR3)*Y/LSPAN  +  LR3 
T3  =  (TS  PAN-TROOT)  *  Y /LSPAN  +  TROOT 
TLE  =  (TSLE-TRLE)  *  Y/LSPAN  +  TRLE 
MZ3  =  (T3-TLE)/(2.  *B3) 

BZ3  =  TLE/2. 

C  CALCULATE  LEADING  SECTION  INERTIA 
C  (RECTANGULAR  PORTION) 

XI3R  =  (1/12.)*B3*(TLE*  *3) 

C  (TRIANGULAR  PORTION) 

XI3T  =  B3*(((MZ3*B3)**3)/6.  +  (2.*BZ3/3.)*(MZ3*B3)**2  + 

&  MZ3#B3*BZ3**2) 

C  TOTAL 

XI3  =  XI3R  +  XI3T 

C  TOTAL  INERTIA  OF  CROSS-SECTION 
XITOT  =  XII  +  XI2  +  XI3 

C  *****  CALCULATE  BENDING  STIFFNESS  ***** 

C  CALCULATE  MOMENT/DYNAMIC  PRESSURE  AT  Y  AND  CHANGE  IN 

MOMENT  M0M2  =  ATOT*(Y-LSPAN)-LRTOT/2.*(Y**2.-LSPAN**2.)- 
&  (LTTOT-LRTOT)*(Y**3.~LSPAN**3.)/(6.*LSPAiN) 

DELMOM  =  MOM2-MOM1 

C  CALCALUTE  THETA/DYNAMIC  PRESSURE  AT  Y  AND  CHANGE  IN  THETA 
RAT  =  MOM2/XITOT 

SUMMOI  =  SUMMOI  +  (MOM2/XITOT)*DY 
THETA2  =  SUMMOI/E 
DELTHT  =  THETA2  -  THETA1 
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C  CALCULATE  DEFLECTION/DYNAMIC  PRESSURE 
ZDEF  =  ZDEF  +  THETA2*DY 

C  CALCULATE  FIN  STIFFNESS 

KFINY  =  AHS(DELMOM  DOITHT) 

MOM1  =  MOM2  THETA1  =  THETA2 

C  *****  WRITE  DATA  TO  FILE  ***** 

C  PRINT  EVERY  10th  DATA  POINT 
IF  (JPRINT  .EQ.  10)  THEN 

C  WRITE  SPRING  RATE  AT  Y  TO  FILE 

WR1TE(4,  100) Y,  KFINY,  MOM2,  THETA2,  ZDEF 
100  FORMAT(2X,  F6.4,  2X,  E15.3,  2X,  FI  1.3,  2X,  E15.3,  2X,  E15.3) 

JPRINT  =  0 
ENDIF 

JPRINT  =  JPRINT+1 
10  CONTINUE 
END 
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DISTRIBUTION  LIST 


copies 


AMSMI-RD  1 

AMSMI-RD-C'S-R  15 

AMSMI-RD-CS-T  1 

AMSMI-GC-1P.  Mr.  Fred  M.  Bush  1 

U.  S.  Army  Materiel  System  Analysis  Activity  1 


ATTN:  AMXSY-MP  (Herbert  Cohen) 
Aberdeen  Proving  Ground,  MD  21005 


IIT  Research  Institute  1 

ATTN:  GACIAC 

10  W.  35th  Street 

Chicago,  IL  60616 

AMSMI-RD-GC,  Mr.  Sliz  1 
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